Quantum Monte Carlo study of the formation of molecular polarizations and 
the antiferroelectric ordering in squaric acid crystals 
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Effects of geometrical frustration and quantum fluctuation are theoretically investigated for the 
proton ordering in a quasi-two-dimensional hydrogen-bonded system, squaric acid crystal. We elu- 
cidate the phase diagram for an effective model, the transverse-field Ising model on a frustrated 
checkerboard lattice, by using quantum Monte Carlo simulation. A crossover to liquidlike paraelec- 
tric state with well-developed molecular polarizations is identified, distinguishably from long-range 
ordering. Emergence of long-range order from the liquidlike state exhibits peculiar aspects originat- 
ing from the lifting of quasi-macroscopic degeneracy, such as colossal enhancement of the transition 
temperature and a vanishingly small anomaly in the specific heat. 

PACS numbers: 77.80.-e,77.84.Fa,75.10. Jm,75.40.Mg 



I. INTRODUCTION 

Proton ordering in hydrogen-bonded systems has long 
been one of the central topics in condensed matter 
physics. Each proton is in a double-minimum potential 
on the hydrogen bond, and spatial correlations among the 
proton configurations strongly affect macroscopic prop- 
erties of hydrogen-bonded crystals. A famous example 
is the "ice-rule" configuration of protons and its rela- 
tion to the residual entropy in water ice [J . Another 
example is the ferroelectricity due to the proton order- 
ing in KH2PO4 Electronic polarization emerging 
from proton ordering has been extensively studied with 
emphasis on both the fundamental physics and the ap- 
plication to electronic devices @, @] ■ 

Here, we focus on one of such hydrogen-bonded ma- 
terials, squaric acid crystal H2C4O4 (H 2 SQ). H 2 SQ is a 
quasi-two-dimensional (2D) molecular solid. In each 2D 
layer, squaric acid molecules form a network of hydrogen 
bonds, as shown in Fig. [lja) Q. A particular configura- 
tion of protons induces a polarization in each molecule, 
and an ordering of the polarizations can lead to ferroelec- 
tricity. In fact, H2SQ exhibits antiferroelectricity below 
T c = 375K, which is driven by the 2D ferroelectric order- 
ing with interlayer antiferroelectric coupling [9| . 

H 2 SQ has two striking aspects. One is the local con- 
straint on the proton positions similar to the ice rule in 
water ice. Each molecule has four hydrogen bonds; two 
out of four protons come close to the C4O4 unit and the 
other two are far. The local constraint alone is not suf- 
ficient to determine a unique ground state and brings 
about a macroscopic degeneracy, as seen in water ice; 
H2SQ has geometrical frustration in nature. 

The other aspect is the effect of quantum tunneling of 
protons. In general, the external pressure increases the 
tunneling rate of protons between two potential minima, 



which reduces local polarizations, and consequently, sup- 
presses the ferroelectricity. Indeed, in H2SQ, T c is sup- 
pressed with increasing pressure. A peculiar intermedi- 
ate state, however, appears before the polarization is lost 
in each molecule: the macroscopic polarization vanishes 
while the polarization in each molecule is retained [Ioj |. 
Consequently, a quantum paraelectric state is realized in 
the low-temperature limit. 

There have been many theoretical studies for the an- 
tiferroelectric transition in H2SQ. The local constraint 
and associated frustration were considered on the basis 
of vertex models or frustrated pseudospin models [TlTtl3 | . 
The cou pling b etween pseudospins and phonons was also 
studied [15l-[l7|. Most of the studies, however, were lim- 



ited at the mean-field level, and the effect of geometrical 
frustration has not been fully clarified yet. In particu- 
lar, quantum fluctuation under the geometrical frustra- 
tion, which is presumably important for understanding 
the quantum paraelectricity under pressure, has not been 
seriously considered so far. 

In the present study, we investigate an effective model 
for H2SQ, a 2D checkerboard-lattice Ising model with 
transverse field which corresponds to the application of 
external pressure. With a sophisticated quantum Monte 
Carlo (QMC) method, we map out the numerically-exact 
phase diagram. We identify a liquidlikc state interven- 
ing between the ferroelectric phase and the paraelectric 
phase. In the intermediate state, molecular polarizations 
are well retained by ice-rule type local correlations, but 
they are globally disordered. The peculiar nature of tran- 
sition from the intermediate state to ferroelectric phase 
is discussed. 

The organization of this paper is as follows. In Sec. [Til 
we introduce models and methods. After introducing the 
pseudo-spin model in Sec. Ill Al we explain the numerical 
method and the definition of the observables in Sec. Ill Bl 



2 




FIG. 1. (Color online) (a) Schematic picture of a layer of 
H2SQ molecules. The dashed lines represent hydrogen bonds 
connecting C4O4 units, and small circles on the bonds show 
protons, (b) Pseudospin representation of the proton dis- 
placement in (a). Ji, J2, and J3 denote the interactions in 
Eq. {TJ. (c) An example of intermediate liquidlike states in 
which every molecule bears a polarization but the system is 
globally disordered, (d) A ferroelectrically-ordered state. In 
(c) and (d), the bold arrows in the center of plaquettes rep- 
resent the molecular polarizations. See the text for details. 

and Sec. Ill CI respectively. The results of calculation are 
presented in Sec. IIIII The temperature dependence of 
observables is given in Sec. IIII Al and the phase diagram 
in Sec. IIII 51 Discussions on our results with the previous 
studies are elaborated in Sec. [TV] Section [V] is devoted 
to summary. 



II. MODEL AND METHOD 

A. Pseudo-spin model 

We here consider a pseudospin model for H 2 SQ fol- 
lowing the previous studies 0, [l4| • In the pseudospin 
model, proton displacements in a plane of the bipartite 
square lattice of H 2 SQ molecules [Fig. [TJa)] are repre- 
sented by the ^-component of pscudospins as of = ±1 
[Fig. [TJb)] [H, Here the appropriate signs are as- 
signed so that protons belonging to A(B) sublattice H2SQ 
molecules correspond to the up (down) spins. The lo- 
cal constraint similar to the ice rule in water ice [l|, Q 
(two out of four protons are close and the other two are 
far) is taken into account by the antiferromagnetic in- 
teractions between nearest neighbors, Ji, and crisscross- 
ing next-nearest neighbors, J2, on the checkerboard lat- 
tice [TH, [HI . When Ji = J2, the model is a 2D variant 
of the spin-ice model [lH, [l9| , in which sixfold degener- 



acy in each plaquctte results in a macroscopic number of 
energetically-degenerate ground states ~ 1.5 N ' 2 . In the 
present case, we take J2 > J\ since two closer protons fa- 
vor an edge of the C4O4 square, not a diagonal fl3l. fl4|. 
J2 larger than J\ partially lifts the degeneracy, but the 
ground-state degeneracy still remains: all configurations 
with different stacking of antiferromagnetic ID diago- 
nal chains, exemplified in Figs. [IJc) andQJd), give the 
same lowest energy. The degeneracy is reduced but still 
quasi-macroscopic, i.e., 4^ . Hence, the present J1-J2 
model does not show any long-range order down to zero 
T [ll|, [l2j , although a finite-T transition was discussed 
in the previous mean- field studies [HI, [l4| . 

To stabilize a ferroelectric ordering [a stripe ordering 
in terms of pseudospins, shown in Fig. Hid)], a further 
degeneracy-lifting perturbation must be included. In the 
present study, we consider an inter-molecular coupling 
originating from the distortion by forming C=C double 
bonds. A C=C double bond induces a trapezoid-typc 
distortion of C4 square, and favors a ferro-type align- 
ment of the molecules since a short C=C bond tends to 
elongate the neighboring parallel C-C bonds in the adja- 
cent molecules. The inter-molecular correlation is incor- 
porated by a third-neighbor ferromagnetic interaction J3 
[Fig. Hfb)]. J3 lifts the remaining degeneracy and selects 
the ferroelectrically-ordered state. 

By summarizing the above argument, our model is 
given by the following Hamiltonian 

H = Jil>K + J 2 ]Ta.f - J 3 E + r E<*) 

where, af is the a component of the Pauli matrix, repre- 
senting the pseudospin operator at site i; J\, J2, J3 > 0, 
and the first, second, and third sums are taken between 
the nearest, second (crisscrossing), third neighbors on the 
checkerboard lattice, as shown in Fig.[ljb). Here, the last 
term with the transverse field V is introduced to repre- 
sent the quantum tunneling of protons. We take J\ = 1 
as the energy unit, and focus on the case with J 2 = 2. 
The results are qualitatively the same for J 2 > J\. We 
analyze the 2D model by changing J3, Y, and T to clarify 
the nature of inplane ferroelectricity in H2SQ; the effect 
of the interplane coupling will be mentioned later. 

B. Quantum Monte Carlo method 

To investigate the thermodynamics of the model 
given by Eq. ([T]), we employ a recently-developed 
continuous-time QMC method with a cluster update in 
the imaginary-time direction (20| . In order to overcome 
the slow relaxation in the present frustrated system, we 
use the replica exchange method [2l| and the loop-flip 
update algorithm [22j ■ In the replica exchange, the repli- 
cas are chosen along the constant Y/T lines since the 
Boltzmann weight strongly depends on the numbers of 
domain walls along the imaginary-time direction which 
are proportional to Y/T. For the loop-flip algorithm, we 
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applied two variants of the original algorithm. One is the 
original algorithm which forms a loop by connecting up 
spins and down spins alternatively, with treating the J\ 
and J2 bonds equivalently. The other one is the diago- 
nal flipping process; a diagonal chain of spins connected 
by J 2 is flipped at once. The system sizes are taken 
from L = 24 to 36 (N = AL 2 ) under periodic boundary 
conditions. Typically, MC measurements are performed 
for 100000 samplings after 30000 initial thermalizations. 
Results are divided into six bins to estimate statistical 
errors by the variance among the bins. 



C. Physical quantities 

To distinguish the long-range order and "ice-rule" type 
local correlations, we calculate the macroscopic polariza- 
tion P and local correlation parameter p. P is calculated 
as 
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P = [S(0, it) 2 + S(tt, 0) 
via the spin structure factor 

S M = iiYl a * a J ex pH k • r *j). 



(2) 



(3) 



which detects the stripe-type ordering in Fig. QJd). The 
critical temperature is determined by the Binder analy- 
sis p3| using the Binder parameter for P, 
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On the other hand, p detects the fourfold-degenerate sta- 
ble configuration in each plaquette by 



P 



(5) 



where the sum p runs over all the crisscrossing plaque- 
ttes, and f(p) is a function giving 1 for the fourfold stable 
states and otherwise —1/3 in pth plaquette: p — > when 
the spins are completely disordered, and p — > 1 when all 
the plaquettes are in the fourfold stable states. Note that 
p is not an order parameter but characterizes a crossover 
associated with the formation of molecular polarizations 
as demonstrated below. We also measure the correspond- 
ing susceptibilities, xp an d Xpi from the fluctuations of 
P and p as 



xp = ^({p 2 )-{p)% 
xp = f «p 2 )-(p> 2 ), 

respectively. The specific heat is calculated by 
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FIG. 2. (Color online) QMC results of (a)(e) P and X p (xp 
is divided by 1000 and 500, respectively), (b)(f) p and X p, an d 
(c) (g) C for the system sizes iV = 4L 2 ranging from L — 24 to 
L = 36. The data are at (a)-(c) J 3 = 0.002 and (e)-(g) J 3 = 
0.004. (d) and (h) shows the Binder parameter in vicinity of 
the T c for J3 = 0.002 and J3 = 0.004, respectively. Insets of 
(c)(g) show the low-T behaviors of C/T. All the results are 
calculated at J2 = 2 along the axis with T/T = tan(7r/6). See 
the text for details. 



III. RESULTS 
A. Locally-correlated liquidlike state 

Figure [U shows QMC results along T/T = tan(7r/6) 
axis at J 3 = 0.002 and 0.004. At the lowest T, the sys- 
tem shows a ferroelectric ordering with fully-saturated 
polarization P, as shown in Figs. [UJa) and[^e). With 
increasing T and T, P steeply decreases and the cor- 
responding susceptibility xp shows a sharp peak which 
grows as N increases, indicating a phase transition into 
a paraelectric state. The critical temperature T c is es- 
timated from the crossing point of the Binder parame- 
ter of P given by Eq. (H|), as shown in Figs. [UJd) and 
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12(h): T c = 0.60(30) for J 3 = 0.002 and T c = 0.76(12) 
for J 3 = 0.004. On the other hand, the local correlation 
parameter p remains to be large even above T c and grad- 
ually decreases for T > 1.0 [Figs. [2(b) and [2(f)]. Corre- 
spondingly, the specific heat C divided by T has a peak, 
as shown in Figs. [2(c) and [2(g). The susceptibility for 
p, Xp, shows a broad peak at a higher T. These indicate 
that the system does not directly enter into a completely 
disordered state at T c but exhibits an intermediate state 
in which molecular polarizations are retained; the system 
shows a crossover to a completely-disordered paraelectric 
state characterized by the peak of C/T or \p a * ^c/t or 
T* p . The intermediate state is a liquidlike paraelectric 
state originating from the ice-rule type local correlations. 

In the ferroelectric phase transition and crossover to 
liquidlike phase, most of the entropy is released at the 
crossover by forming the ice-rule type manifold. As 
shown in Figs. [2(c) and [2(g), C/T sharply decreases 
with the saturation of p, and the peak associated with 
the phase transition at a lower T is very small. In fact, 
as shown in the insets, the intensity of the small peak 
decreases as N increases, while the peak position ap- 
proaches T c . This suggests that the entropy associated 
with the phase transition becomes vanishingly small in 
the thermodynamic limit. The peculiar behavior is un- 
derstood by considering the quasi-macroscopic degener- 
acy in the ice-rule type manifold where the remaining 
entropy is in the order of y/~N not N [l2T |. 

The intermediate liquidlike state is widely observed 
while changing T/T. Figure [3] shows P, p, and their 
susceptibilities along the various T/T axes. As indicated 
by a decrease of p and a broad peak of \ p , the crossover 
temperature T* largely decreases with increasing T, while 
the data of P and \p show that T c does not decreases 
so rapidly. There, however, is always a window of the 
intermediate state with well-developed local correlations 
and suppressed global order. 
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FIG. 3. (Color online) T dependences of (a) P and p, and 
(b) x^/2000 and Xpi measured along various T/T. T/T is 
parametrized by 9 (degree) with T/T — tan 9. All the results 
are calculated at J2 = 2 and Jz = 0.002 for N = 4 x 36 2 sites. 



B. Phase diagrams 

The phase diagrams for T, T, and J3 are summarized 
in Fig. [4] (T C ,T C ) are estimated by the Binder analy- 
sis of P, and (T*,r*) are identified by a peak of \ p or 
C/T. In the case of J3 = 0, the system exhibits only the 
crossover into the liquidlike state with quasi-macroscopic 
degeneracy [the region below T* ; see Figs. 2(b) and (f)], 
and remains paraelectric down to the lowest T calculated. 
One might expect a quantum order by disorder as pre- 
dicted for the case of J\ = J2 but it will be limited 
to very low T region, if any, and is out of the scope of 
the present study. With J3 switched on, the degeneracy 
is lifted and the ferroelecrically-ordered phase emerges 
inside the liquidlike state. The ordered state rapidly ex- 
tends with increasing J3, but a sequential change of the 
three different regimes is clearly observed in the region 
where J3 is sufficiently small. With further increasing 
J3, T c exceeds T* and the intermediate liquidlike state is 



taken over by the long-range ordered state. 

A remarkable point of the phase diagram is the rapid 
growth of T c with J3, which is more than 100 times larger 
than J 3 . This colossal enhancement of T c is also as- 
cribed to the peculiar nature of the intermediate state 
with quasi-macroscopic degeneracy; the spatial correla- 
tions are strongly enhanced along the diagonal chains. 
This strong correlation reinforces the effect of J3 as J2 
and J3 are not frustrated, and the system becomes ex- 
tremely sensitive to the degeneracy-lifting perturbation 
J3. Another point is the fate of the intermediate state 
at low temperatures. Since the QMC simulation be- 
comes harder for larger T/T, it is difficult to conclude 
whether the intermediate state remains at low tempera- 
tures. Nonetheless, we expect a finite window down to 
low T for small J3 from the systematic change of the 
phase diagram shown in Fig. 4. 
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IV. DISCUSSIONS 
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Comparing with experiments, we successfully identify 
the liquidlike state with well-developed molecular polar- 
izations, which might account for the peculiar intermedi- 
ate state in experiments [l(| • Our results suggest that the 
intriguing physics related with the ice-rule type degener- 
acy is involved within the 2D layers of H2SQ crystal. We 
note, however, that the qualitative shape of the phase 
diagram appears to be different; in experiments, both T c 
and T* decrease almost linearly in applied pressure and 
their difference is almost independent of pressure. This 
can be ascribed to the parametrization of the realistic 
situation; i.e., how the model parameters change under 
pressure. The first-principle calculations may help fur- 
ther quantitative studies [25[ . Another possible origin of 
the discrepancy is an ambiguity in assigning T*; in gen- 
eral, crossover boundary depends on how to define or de- 
tect it. With regard to the nature of the phase transition, 
our results indicate that it is second order and the asso- 
ciated anomaly of the specific heat is vanishingly small. 
These apparent contradictions with experiments [26M28| 
might be reconciled by considering the intcrlaycr cou- 
pling or more complicated couplings to lattice distortions, 
which arc neglected in our model 

Comparing with the previous theoretical researches, 
our results are obtained by seriously including both ge- 
ometrical frustration and quantum fluctuation, which 
were not fully taken into account in the previous studies. 
For instance, in the absence of J3, our result shows no 
phase transition as expected in the frustrated situation, 
in contrast to the previous research in which a finite- 
temperature phase transition was predicted because of 
the mean-field type treatment. Furthermore, our result 
clearly indicated existence of the liquidlike state as well as 
several significant consequences of the local correlation on 
the thermodynamical properties. In particular, the ab- 
sence of anomaly in the specific heat at T c and the strong 
enhancement of T c by J3 are revealed for the first time 
by our calculations. On the other hand, as mentioned 
above, our result also showed continuous transition as 
in many previous theoretical studies using pseudospin or 
an equivalent approach. This suggests that an extension 
of the model is necessary to account for the first-order 
transition in experiments. 



SUMMARY 



FIG. 4. (Color online) Phase diagrams obtained for vary- 
ing J3. Each diagram for (a)-(d) corresponds to J3 = 
0, 0.002, 0.004, 0.008, respectively. The solid circles shows T c 
obtained by Binder analysis and triangular (square) points are 
T* evaluated from C/T (Xp)- The dotted line shows phase 
boundary of ordered state and paraelectric state, and the in- 
termediate state is assigned as the region between T£ p and 
T c . The gradation shows intensity of p. 



To summarize, we have investigated the geometrically- 
frustrated transverse-field Ising model as an effective 
model for squaric acid crystals. By unbiased QMC sim- 
ulations, we have identified an intermediate liquidlike 
state between the ferroelectrically-ordered state and the 
completely-disordered paraelectric state, in which molec- 
ular polarizations are well preserved but they are glob- 
ally disordered due to the frustration. Furthermore, we 
found out that the emergence of locally correlated state 
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significantly affects the thermodynamic behavior of this 
system. In particular, we unveiled the vanishingly-small 
anomaly in the specific heat at the transition point and 
the colossal enhancement of T c by the degeneracy-lifting 
perturbation J3. Emergence of such state and the re- 
markable effects have never been reported in the previ- 
ous theories. The liquidlike state accounts for a peculiar 
intermediate paraelectric state observed under external 
pressure in the squaric acid crystal. 

Although our results qualitatively reproduce the pe- 
culiar phase diagram of the squaric acid under pressure, 
the quantitative changes of the critical temperature and 
the crossover temperature are different from those in ex- 
periments. Further quantitative researches are neces- 
sary. For example, experimentally, detailed analysis of 
the structure under pressure will be quite important for 
more quantitative comparison between experiment and 
theory. Moreover, the first-principle calculation will help 
to identify the model parameters and their changes under 



pressure more precisely. On the other hand, our result 
exhibits second order phase transition, which is in con- 
trast to the first order transition observed experimentally. 
This might require further extension of the model, for 
example, by including the intcrlayer coupling and com- 
plicated couplings to lattice distortions. Such extension 
is left for the future work. 
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